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Abstract 

We investigate the jamming transition in packings of emulsions and granular materials via molec- 
ular dynamics simulations. The emulsion model is composed of frictionless droplets interacting via 
nonlinear normal forces obtained using experimental data acquired by confocal microscopy of com- 
pressed emulsions systems. Granular materials are modeled by Hertz-Mindlin deformable spherical 
grains with Coulomb friction. In both cases, we find power-law scaling for the vanishing of pressure 
and excess number of contacts as the system approaches the jamming transition from high volume 
fractions. We find that the construction history parametrized by the compression rate during the 
preparation protocol has a strong effect on the micromechanical properties of granular materials but 
not on emulsions. This leads the granular system to jam at different volume fractions depending on 
the histories. Isostaticity is found in the packings close to the jamming transition in emulsions and 
in granular materials at slow compression rates and infinite friction. Heterogeneity of interparticle 
forces increases as the packings approach the jamming transition which is demonstrated by the 
exponential tail in force distributions and the small values of the participation number measuring 
spatial localization of the forces. However, no signatures of the jamming transition are observed 
in structural properties, like the radial distribution functions and the distributions of contacts. 



I. INTRODUCTION 



A variety of systems such as granular materials, compressed emulsions, or molecular 
glasses, exhibit a non-equilibrium transition from a fluid-like to a solid-like state. When 
the constituent particles of such systems are so crowded that they are in close contact with 
one another, the whole system experiences a sudden dynamical arrest, which is referred 
to as the glass or jamming transition . A jammed system is a many-body system 

blocked in a configuration far from equilibrium, from which it takes too long a time to relax 
in the laboratory time scale. Jammed systems are usually random packings of particles, 
which have been the starting point for many studies of liquids and glasses j&U, ^,0, Z|- A 
;heoretical approached to treat these systems has been developed by Edwards and coworkers 
^, 0] who have proposed a statistical mechanics description of such jammed system, in 
which thermodynamic quantities are computed as flat average over jammed configurations 

In a recent paper [13] we found that the jamming transition in packings of granular 
materials and emulsions can be viewed as a phase transition at the critical concentration of 
random close packing (RCP) for frictionless emulsions or random loose packings (RLP) for 
frictional grains. The transition is characterized by a power behavior and the corresponding 
set of exponents describing the stress, coordination number and elastic moduli as the system 
approaches the jam ming transition. This idea has also been investigated for packings of 
frictionless particles Ilfij] and colloids Here we elaborate further on the results 

of |3| and present our detailed computer simulation studies about many aspects of the 
jamming transitions in these systems. 

We consider packings of spherical deformable particles interacting via friction and friction- 
less forces. Granular packings are characterized by elastic Hertz-Mindlin forces 18. Iroj] 
with Coulomb friction. Concentrated emulsion systems are modeled as a collection of fric- 
tionless deformable droplets interacting via a nonlinear normal force law which is obtained 

n 

from the analysis of experimental data acquired using confocal microscopy in j2(lj . 

We first introduce a protocol to generate a series of packings near the jamming transition 
and beyond for both emulsions and granular materials. Then we study how the construction 
history of the packings affects the jamming transition. The existence of frictional tangential 
forces is one of many generic properties of granular materials. In the elastic Hertz-Mindlin 
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model the tangential force is proportional to the tangential displacement, and it is trun- 
cated by a Coulomb cut-off. Both the elastic tangential force and Coulomb friction are 
path-dependent. The path-dependence in the microscopic interactions leads to the path- 
dependence of the granular packings at the macroscopic level, an effect that has been shown 
in several classical experiments P-kI I2I | . In our simulations we find that the construction 
history affects the properties at the jamming transition in the granular packings through 
their direct impact on the tangential forces, while in the case of emulsions we do not find 
any path-dependence. 

Next, we fully characterize the evolution of the micromechanical and microstructural 
properties, such as the stress, the coordination number and isostaticity, the radial distribu- 
tion function, the probability distributions of normal and tangential forces and interparticle 
contacts, the existence of force chains, the participation number and plasticity index, as the 
system approaches the jamming transition. We found that the heterogeneity of interparticle 
forces increases as the packings approach the jamming transition which is demonstrated by 
the exponential tail in force distributions and the small values of the participation number. 
However, no clear signatures of the jamming transition are observed in structural properties, 
like the radial distribution functions and the distributions of contacts. 

The paper is organized as follows: in Section II we discuss the microscopic force models 
used in our simulations of granular materials and emulsions. In Section III we give the 
details of the numerical protocol to generate packings near the jamming transition and 
beyond and the various measurements to be performed. The main results of this work are 
given in Section IV and we summarize our results in Section V. 

II. MICROSCOPIC MODEL 

In this section we discuss the microscopic models of inter-particle forces for granular 
materials and emulsions. We first briefly review the standard Hertz-Mindlin model for 
granular materials; more details can be found in ^ 

" Princen model" |23[ for the inter-particle forces in emulsions and its nonlinear modification 
according to experimental data on deformation in concentrated emulsion systems 



IS , l23|. Then we describe the linear 
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A. Contact mechanics in granular materials 



We describe the microscopic interaction between grains by the nonlinear Hertz-Mindlin 
normal and tangential forces. The normal force between two contacting grains at position 
~x\ and ~x 2 with uncompressed radii R\ and R 2 is 1?], 0| 

Fn = IknR^e' 2 , (1) 

where R is the geometric mean of R\ and R 2 : R — 2R\R 2 / [R\ + R2), and £ is the normal 
overlap: £ = (1/2) [Ri + R 2 — \ x 1 — x 2 \] ■ The normal force acts only in compression, i.e. 
F n = when £ < 0. The effective stiffness k n = AG/ (1 — v) is defined in terms of the shear 
modulus of the grains G and the Poisson ratio v of the material the grains are made. From 
Eq. we can see that the normal forces are completely determined by the geometrical 
configuration of the packing and have nothing to do with the construction history. 

The tangential contact force was first calculated by Mindlin [18] for grains under oblique 
loading. In his model, the increment in tangential force is 

AF t = k t (R0 1/2 As (2) 

where k t = 8G/(2 — u), and the variable s is defined such that the relative tangential 
displacement between the two grain centers is 2s. Therefore the tangential force is obtained 
by integrating Eq. (J2J) over the path taken by the particles in contact with the initial 
condition, F n = 0, F t = at £ = 0, s = 0, yielding: 

k t (R0 1/2 ds . (3) 

path 

As the tangential displacement increases, the elastic tangential force F t reaches its limiting 
value given by a Coulomb cut-off for granular materials, which is 

F t = nF n . (4) 

The Coulomb cut-off adds a second source of path-dependence to the problem. 

From the discussion above, one can see that a straightforward way to affect the tangential 
force is to change the frictional coefficient fi, which controls the maximum value of the ratio 
between tangential and normal force. But the elastic tangential force, Eq. (j2J), is even 
more intimately related to the tangential displacement, which is mainly determined by the 
construction history. 



B. Force law for emulsion 



Emulsions are a class of material which is both industrially important and exhibits very 
interesting physics [23]. They belong to a wider material class of colloids in that they consist 
of two immiscible phases one of which is dispersed into the other, the continuous phase. Both 
of the phases are liquids and their interface is stabilised by the presence of surface-active 
species. Emulsions are composed of droplets of a liquid (for instance silicone oil) stabilized 
by a surfactant (like sodium dodecylsulphate ) in a continuum phase (such as a water and 
glycerol mixture) j^j]. Being composed of only liquids, emulsion droplets interact with each 
other only via normal forces with no solid friction between them. The determination of an 
accurate force model for the compression of two droplets is by no means trivial, but it can be 
simplified in certain limits. For small deformation with respect to the droplets surface area 
the Laplace pressure remains unchanged and all energy of the applied stress is presumed to 
be restored in the deformation of the surface. Then the normal repulsive force between two 
spherical droplets in contacts with uncompressed radii R\ and R2 can be calculated at the 
microscopic level as: 

/ = J^Amea (5) 

This is known as the Princen model j^J, where a is the interfacial tension of droplets, / is 
the normal force between two droplets, R is the geometric mean of R\ and R 2 , and A mea is 
the measured contact area between the droplets. As in the granular materials, the normal 
force acts only in compression. 

In j^j] a linear model was used to relate the area of deformation with the overlap between 
the spheres £, resulting in a linear spring model for the force-law between the droplets. More 
detailed calculations j^J and numerical simulations j^j show that the interdroplet forces in 
emulsions are better represented by a nonlinear spring / oc £ a with the exponent a between 
1 and 1.5 and, more importantly, depending on the number of contact forces acting on the 
droplet. 

Recently, Brujic et al. j.2lj used confocal microscopy to study a compressed polydisperse 
emulsion. This system consists of a dense packing of emulsion oil droplets, with a suffi- 
ciently elastic surfactant stabilising layer to mimic solid particle behaviour, suspended in a 
continuous phase fluid. By refractive index matching of the two phases they obtained a 3D 
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image of the packing structure by using confocal microscopy for the imaging of the droplet 
packings at varying external pressures, i.e. volume fractions. The key feature of this optical 
microscopy technique is that only light from the focal plane is detected. Thus 3D images of 
translucent samples can be acquired by moving the sample through the focal plane of the 
objective and acquiring a sequence of 2D images. One 2D image obtained experimentally is 
shown in Fig. QJi. 

The emulsion system, stable to coalescence and Ostwald ripening, consisted of Silicone 
oil droplets (rj = lOcS) in a refractive index matching solution of water (w t = 51%) and 
glycerol (w t = 49%), stabilised by 20mM sodium dodecylsulphate (SDS) upon emulsification 
and later diluted to below the critical micellar concentration (CMC= 13mM) to ensure a 
repulsive interdroplet potential. The droplet phase is ffuorescently dyed using Nile Red, 
prior to emulsification. The control of the particle size distribution, prior to imaging, is 
achieved by applying very high shear rates to the sample, inducing droplet break-up down 
to a radius mean size of 3.4/xm. This system is a modification of the emulsion reported by 
Mason et al. J3] to produce a transparent sample suitable for confocal microscopy. 

In these images, the area of contacts, the droplets and the aqueous background differ by 
the darkness in an 8-bit gray-scale image such as the one displayed in Fig. which has an 
average darkness 210, 90, and 30, respectively. The result of the image analysis carried out 
by Brujic et al. is a set of contact areas, A mea , along with the undeformed radii (R\, R^) and 
positions (xi, £2) of two droplets giving rise to each contact. Using the obtained information 
of Ri and ~r \ and ignoring the droplet deformations, we can reconstruct the images as shown 
in Fig. [TJd and calculate a geometric overlapping area 

A ca i = TLRf , (6) 

where £ is the overlapping. For all the 1439 contacts obtained in the experiments, A ca i is 
plotted against A mea in the inset of Fig. EJ from which we can see that A ca[ is generally 
different from A mea , indicating that there exist deviations from the linear force law Eq. (|5|). 
This fact provides a direct measure of the effects of anharmonicity of interaction between 
the droplets surfaces. 

In order to extract the nonlinear dependence of the force law between the force and the 
deformation characterized by £, we plot the calculated area A ca \ from the reconstructed 
image versus the real area of deformation A mea . In order to achieve this, we need to improve 
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(a) (b) Cc) 



FIG. 1: Image slice of emulsion packing from confocal microscope (b); reconstructed image before 
(a) and after Monte Carlo annealing (c). 

the estimation of the droplet centers and radii by minimizing the difference between the 
reconstructed image and the original one. In this minimization, i?, and Xj are the changing 
variables and we use a Monte Carlo annealing method to find the optimum set of parameters. 
The corresponding image after the minimization is shown in Fig. which is closer to the 
original image than the reconstructed image shown in Fig. as indicated by rectangular 
region in the figure. Finally, we choose only those contacts for which Ri and Xi change less 
than 5% before and after the minimization, and plot A ca i versus A mea in the main panel of 
Fig. |21 The plot shows a clear trend which can be fitted as 

A cal = 2.67 (A mea f 78 . (7) 

Combining Eqs. (jHJ), © and (JZJ), we obtain the nonlinear force-law for emulsion droplets: 

/ = 1.23 a R 028 e 28 - (8) 

In our simulations, Eq. (jHJ, which takes in account partially the effects of the anharmonicity 
of interaction between the droplets, is used. We also use a linear spring force law / = 47ro"£, 
which completely neglects the anharmonicity of interaction. These two forces give similar 
results in the quantities we measured. In this paper, we only present the results using Eq. 



7 




FIG. 2: Acai is plotted against A mea after Monte Carlo annealing; the red line is a fit of . In the 
inset, A ca i is plotted against A mea for all forces before Monte Carlo annealing. 



III. SIMULATIONS, PROTOCOL AND MEASUREMENTS 

In this section, we first describe the method of molecular dynamics simulations. Then 
we show the protocol to generate packings approaching the jamming transition, and in the 
third part we define various quantities we use to characterize the jamming transition. 



A. Molecular dynamics simulation 

We follow the Discrete Element Method (DEM) or Molecular Dynamics developed by 
Cundall and Strack [^J and solve the Newton equations for an assemble of particles with 
zero-gravity, which interact via the microscopic models discussed in Section II. DEM em- 
ploys a time-stepping, finite-difference approach to solve the Newton equations of motion 
simultaneously for every particle in the system: 

^ d 2 X 



ni- 



di 2 



(9) 
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M = I— (10) 

where F and M are the total force and torque acting on a given particle, m and / are the 
mass and moment of inertia, X and 6 are the position and angle of the particles respectively. 
The numerical solution of Eq. Q and Eq. fllOl) are obtained by integration, assuming 
constant velocities and accelerations for a given time step. We choose the time step to be 
a fraction of the time that it takes for sound wave to propagate on one grain (or droplet). 
A global damping proportional to the translational velocity was used to dissipate energy in 
both cases. The damping simulates the drag force of the droplets in the suspending medium. 
More details about simulations can be found in 

In the simulations of emulsions, the system consists of 10000 droplets, which are all 2 
fxm in diameter. Droplets interact via the normal force given by Eq. (jHJ) with the interfacial 
tension a = 9.8 x 10~ 3 N/m and the density p = 10 3 kg/m 3 . The granular system is composed 
of 10000 glass beads of equal size (radius=0.1mm) interacting via the Hertz-Mindlin forces, 
Eqs. (0) and (jHJ) and the Coulomb force Eq. (J3J). The microscopic parameters defining the 
glass beads are the shear modulus: G = 29Gpa, the Poisson's ration: v = 0.2, the friction 
coefficients \x = 0.3 and the density: p = 2 x 10 3 kg/m 3 . 



B. Numerical protocol 

We start our simulations from a set of non-overlapping particles located at random po- 
sitions in a periodically repeated cubic box with an initial volume fraction around <fi ~ 0.2. 
The box is compressed isotropically by constant compression rate 7 until a given volume 
fraction, 0, is reached. Then the compression is turned off and the system is allowed to 
relax with constant volume until it reaches a stable state, which means that the pressure of 
system remains unchanged over a period of time (usually 5 x 10 5 MD steps). This protocol 
can generate packings with different volume fractions. From simulations, we find that there 
is a critical volume fraction C , below which a jammed packing with nonzero pressure can 
not be obtained. This fact is illustrated in Fig. |3J where the time evolution of the pressure 
in this compress-and-relax process for two packings are shown. The packing with 0i > C 
stabilizes at none-zero pressures but the pressure decreases very fast to zero for 02 < <Pc 
even though 0i and 02 are different only by 1.6 x 10~ 5 . 
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FIG. 3: Time evolutions of pressure for two packings in simulation of granular material. Line 
represents a packing with <f>\ > (f) c and circles represent a packing with <f>2 < 4>c where = 
02 + 1.6 x 1(T 5 . 



Several other ways of preparing static packings of particles exist in the literature. Conju- 
gate gradient methods were used to study dense random packing of frictionless particles by 
O'Hern et al. A pouring-ball method, which mimics the process of pouring balls into 

a container under gravity, was used by Silbert et al. 0, 3^. A servo mechanisms which 
adjusts the strain in the system to equilibrate the packing at a given stress was used in 19] . 
One of the advantages of the present protocol is that it allows us to generate packings with 
different construction histories by using different compression rates. Thus we can study the 
path-dependency systematically. 



C. Computation of characteristic quantities 

After generating the packings, we compute the following quantities to characterize their 
micromechanical and structural properties: 
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1. Pressure and coordination number 



The macroscopic stress tensor for point contacts in a volume V is given by 

k 

where V is the volume of the system, R k is the vector joining the center of two particles 
in contact k, nt = R k /R k and F k is the total force in contact k. The pressure a is the 
average of three diagonal elements of the stress tensor, i.e. a = (an + a 2 2 + cr 33)/3. The 
coordination number, Z, is the average of contacts per particle in contact network: Z = 
where M is the total contacts and N is the number of total particle in contact network. Due 
to the non-gravity environment in our simulation, some floaters, which have zero contacts 
and don't participate in the contact network, exist. We exclude them when calculating the 
coordination number. We note the floaters were also reported in Ref. 



2. Force distributions, force chains and participation number 

In both emulsion and granular packing, interparticle forces are highly inhomogeneous. 
This is quantified by the probability distribution of the normal and tangential forces calcu- 
lated for emulsions and grains. Moreover, photoelastic visualization experiments j^l and 



simulations 



13. 



38, 



39| have shown that the contact forces are strongly localized along "force 
chains" which carry most of the applied stress. To quantify the degree of force localization, 
we consider the participation number V j^J: 

r=(M&) • (ii) 

Here M is the number of total contacts, and 

M 

Qi = /i/E/i; 
i=i 

where fj is the total force at contact j. From the definition, T = 1 indicates qi = 1/M, for 
all qi and a state with a spatially homogeneous force distribution. On the other hand, in 
the limit of complete localization, T w 1/M, which is essentially zero in our simulation. 
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3. Plasticity Index 

In order to quantitatively characterize the tangential forces in the granular packings, we 



measure the plasticity index following 



3S 



E ^ < 12 > 

which has value between and 1. The maximum means the contact is "plastic" and the 
tangential force is the Coulomb friction; otherwise the contact is elastic and the tangential 
force is the Mindlin elastic force. 



IV. RESULTS 

A. Power-law scaling near jamming. Effect of construction history and Isostaticity 

Using the protocol discussed above, we generate series of packings of emulsions and 
granular materials approaching the jamming transition. In order to study the effect of 
constructing history, different compression rates are used to generate the packings. For the 
case of granular materials, we use four compression rates: 2 x 10 4 , 2 x 10 3 , 2 x 10 2 and 
2 x 10 1 m/s, while for emulsions we use two: 1.5 x 10~ 2 and 1.5 x 10~ 3 m/s. In Fig. |3]we 
plot the pressure, cr, and coordination number, Z, of the resulting packings as a function 
of volume fraction, <ft, where different symbols correspond to different compression rates. 
Quantitative analysis of the data in Fig. 0] shows that in both case the pressure a and extra 
contacts Z — Z c vanishes as power-laws as cj) approaches the jamming transition at <p c as 
shown in [is] ]: 

(13) 

and 

z - z c ~ (</> - 4> e y. (i4) 

The power-law fits are shown by lines in Fig. H] and the extracted fitting parameters 
including the critical volume fraction <p c , the critical (minimal) coordination number Z c and 
two power indexes a and (3 are listed in Table I. 

As shown in Fig. the jamming transition for emulsions occurs at <ft c = 0.645 and Z c = 6 
regardless of the compression rate, i.e. the construction history. We identify <p c = 0.645 as 
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FIG. 4: Power-law scalings in the packings of granular materials and emulsions approaching the 
jamming transition. The symbols are the data from the simulations and different symbol corre- 
sponds to different compression rate. The lines are power-law fits, (a) a vs <j> in emulsions, (b) Z 
vs 4> i n emulsions, (c) a vs (j) in granular materials and (d) Z vs <j> in granular materials. 



random close packing (RCP) volume fraction, which value has been reported many times 



before 



14J, |32| . The minimal coordination number Z„ = 6 indicates that the packings 
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14, 



31 



near the jamming transition in emulsions are isostatic 

A packing is isostatic (2] when the number of contact forces equal the number of force 
balance equations. In a 3D packing of perfectly smooth (frictionless) particles there are 
NZ/2 unknown normal forces and 3N force balance equations. This gives rise to a minimal 
coordination number needed for mechanical stability as Z c = 6 in 3D (the isostatic limit). 
In the case of packings of perfectly rough particles which is realized by frictional particles 
with infinite friction /i — > oo — notice that in this case there are still tangential forces given 
by the Mindlin elastic component Eq. (j2J) — besides NZ/2 unknown normal forces and 3N 
force balance equations, there are NZ unknown tangential forces and 3N torque balance 
equations. Thus the coordination number in the isostatic limit is Z c = 4 for frictional 
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System 


Com. rate (m/s) 


0C 




a 


P 


Emulsion, fx = 


1.5 x 10" 2 (10~ 3 ) 


0.645 


6.01 


1.25 


0.51 


Granular, /j = 0.3 


2 x 10 4 


0.604 


4.52 


1.52 


0.46 


Granular, fi = 0.3 


2 x 10 3 


0.590 


4.53 


1.46 


0.45 


Granular, = 0.3 


2 x 10 2 


0.581 


4.53 


1.48 


0.46 


Granular, fi = 0.3 


2 x 10 1 


0.576 


4.54 


1.52 


0.47 


Granular, /i = oo 


2 x 10 2 


0.571 


3.98 


1.66 


0.45 



TABLE I: Parameters in the power-law fittings near the jamming transition Eqs. (|13j) and 1)14(1 
for emulsions and granular materials. 



isostaticity is the foundation of recent 
28, 2*1 Q- 



packings in 3D. In such isostatic packings, there is possibly a unique solution for the forces 
between particles for a given geometrical configuration, because the number of equations 
equals to the number of unknowns. The existence o 
theories of stress propagation in granular materials 

Thus, in our simulation the isostatic limit in the frictionless case is approached as a — > 
or — > (j) c . In other words, the isostaticity appears in the limit of rigid balls, or when the 
rigidity of the particles goes to infinity. This corresponds to the so-called marginal rigidity 
state |30j ]. 

Quite different from the results is emulsions, <p c in frictional granular materials with finite 
/j = 0.3 varies depending on the compression rates from 0.576 (slowest rate) to 0.604 (fastest 
rate) as seen in Fig. 0] The critical coordination numbers in Table I for granular materials 
are around Z c w 4.5. This value is above Z c = 4, the value required by isostaticity. However, 
it does not mean that these frictional granular packings with finite \i = 0.3 are not isostatic 
since the prediction Z c = 4 is strictly valid for packings with fi — > oo. We will elaborate 
more on this result below and discuss in more detail this limit. 

Next we try to understand the values of the power indexes. In the granular systems, due 
to the Coulomb cut-off (with fi = 0.3), the intergranular normal forces are always larger 
than the tangential forces. The main contribution to the pressure is from the normal forces 
and not the tangential ones. Therefore, for Hertzian balls we expect the following scaling: 

a „ Fn „ £3/2 „ (0 _ 0c)3 / 2) 



14 



which implies that the exponent a in Eq. (|13p should be equal to the 3/2 exponent of 
the normal force law in Eq. (0). The values of a obtained in our simulations for different 
compression rates shown in Table I are all around 3/2, supporting the above argument. This 
argument should also hold in the case of emulsions since only normal forces act between the 
droplets. In fact, we found a = 1.25 for emulsions independent of the compression rate 
which is close to the exponent 1.28 in the force law Eq. (JHJ). From Table I we find that the 
exponent (3 in Eq. (JTljl are very close to 0.5 for both emulsions and grains, independent of 
the force law and also of the compression rate. This numerical value and the independence 
of the force law were also reported by O'Hern et al. ^j]. These authors provided a possible 
way to understand f3 = 0.5. In [32j it was shown that the radial distribution function, g(r), 
in packings near the jamming transition displays a power-law behavior near r = D, where 
D is the diameter of spheres: 

g{r) oc (1 - r/D)-°- 5 . (15) 

If one assumes an affine deformation upon compression, then one consequence of Eq. f|15|) is 
that the coordination number should increase with the power (3 = 0.5. However, as shown 
below, we do not observe the power-law region in q(r) from our data. Therefore, to our 
opinion, the origin of (3 — 0.5 is still not clear. In [l9| we provide more detailed analysis of 
this argument assuming a narrow peak of the pair distribution function at r = D. However, 
the predicted dependence of the coordination number with pressure still disagree with the 
(3 = 0.5 result. 

Here we elaborate further on the topic and present an analogous derivation of the expo- 
nent (3. Let us assume that, in the limit of zero pressure, there is a probability distribution 
P{h) of gap sizes, h, between each ball and its neighbors: 

P{h) = Z c 5(h) +d 1 hr*, (16) 

which is consistent with a power-law in g{r) as given by Eq. (|T5|) . In ^| we have consid- 
ered a delta function followed by a Taylor expansion around h = instead of the singular 
behavior of Eq. (jlfij) . Following the derivations in Q| which employ an effective- medium 
approximation which assumes an affine motion of the grains with the external perturba- 
tion (being compression of shear) we arrive at the pressure expressed in terms of the static 
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FIG. 5: Distributions of S parameter. The measurements are carried on 4 packings around 650 
KPa generated by different compression rates. 

compressive strain, e < 0, as 

a='^{Z c {-e)^ + d 1 {-e)^-% (17) 

D7T 

and the coordination number becomes 

Z = Z c + d 3 a^ 1 ~^. (18) 

This last result is consistent with Eqs. (fT3|) and (|14jl with (3 /at = (2/3)(l — £). Thus a value 
of £ = 0.5 would fit all our data. However, we could not find evidence of Eq. ()16j) in our 
simulations. 

As we discussed in Sec. II, the tangential forces in the Hertz-Mindlm model are pro- 
portional to the normal displacements and therefore depend on the history of interaction. 
The history-dependence shown in Fig. and0Jl for the granular packings is a result of this 
microscopic path-dependence. In our simulations, a lower compression rate allows the grains 
to have enough time to relax and slide against each other during the preparation protocol. 
This allows for larger tangential displacement between the grains. On the other hand, a 
higher compression rate "freezes" the system so quickly that the grains have no time to 
slide. Thus we expect that the tangential force will be smaller in the packing generated by 
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1 2 3 4 5 6 7 

f t /<f> 



FIG. 6: Probability distribution of the tangential forces in packings generated by different com- 
pression rates. The X and Y axes are scaled by the averaged force and the probability of averaged 
force respectively. This scaling convention will be used in all the probability distribution plots of 
forces throughout the paper. 

fast compression rates. In order to quantify this idea, we measured the plasticity index Eq. 
(J12j) E = -0r- for four packings around 650 kPa generated by different compression rates. 
The probability distributions -P(E), shown in Fig. confirm the above argument. In the 
packing generated by 7 = 2 x 10 1 m/s, almost 20% of the tangential forces are truncated by 
the Coulomb cut-off as evidenced by the fact that P(E) has a maximum at E = 1. On the 
other hand, for the fastest compression rate, 7 = 2 x 10 4 m/s, the peak of the distribution 
is around E ~ 0.4 and P(E = 1) ~ implying that there are almost no plastic contacts. 
This indicates that the grains in the packing generated by lower compression rates feel more 
"sticky" to each other, therefore the packings can sustain non-zero pressure at lower volume 
fraction. This explains why <p c decreases as 7 decreases in granular materials. In Fig. |U1 we 
plot the probability distribution functions of the tangential forces in the packings generated 
by different compression rates. The distributions are different in the small forces regions. 
The packing generated by 7 = 2 x 10 4 m/s contains more small forces, which are less than 
the average forces, than the packing generated by 7 = 2 x 10 1 m/s. 

We have shown that the critical contact number Z c ~ 4.5 in the packings of granular 
materials with /i = 0.3 is larger than the prediction of isostatic considerations Z c = 4. We 
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FIG. 7: a VS Z in the packings of granular material in the infinite friction limit. 



first noticed that Z c = 4 is strictly valid when /x — * oo. Moreover, we showed above that 
a fast compression rate inhibits the developments of tangential forces and therefore may 
also contribute to the breakdown of isostaticity. Therefore, we now investigate the slow 
compression rate limit of a granular packing in the /i — > oo limit. We carried out a set 
of simulations of grains interacting via Hertz forces Eq. (0) and tangential Mindlin forces 
Eq. (J2J but without the Coulomb cut-off (/i = oo) for a compression rate 7 = 2 x 10 2 m/s, 
slow enough to allow sufficient time for grains to rearrange. The result of the simulations is 
shown in Fig. [7| The pressure versus volume fraction can be fitted by 

a ~ (Z - 3.98) 3 69 , (19) 

with a minimal coordination number Z c = 3.98 indicating the recovery of isostaticity. We 
find 4> c = 0.571 and exponents a = 1.66 and (3 = 0.45. We also find that Z c is larger than 4 
if a compression rate larger than 7 = 2 x 10 2 m/s is used. This may also suggests that the 
lack of isostaticity found in (3^ for frictional packings might be due to a fast preparation 
protocol inherent to the pouring-ball method. 

Finally, we note that, besides S and the probability distribution function of tangential 
forces, we also measure the normal force distribution, coordination number distribution and 
radial distribution function for the four packings around 650 kPa at different compression 
rates and found no clear signature of history dependency in these measurements. 
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FIG. 8: Micromechanical properties at different pressures in the packings of granular materials. 
Measurements are carried on the packings generated by 7 = 2 x 10 4 m/s. (a) Distribution of the 
normal forces, (b) Distribution of the tangential forces, (c) T vs pressure and (d) Distribution of 
S. 



B. Micromechanical and microstructural properties 



In order to study the micromechanical and microstructural properties near the jamming 
transition we carry out the measurements described in Section II. B. For the micromechanical 



properties, we measure the normal and tangential force distributions [20 



36 



37 



and the 



T [13( and £ parameters |32( . For the microstructural properties we calculate the distribution 
of contacts and the radial distribution function. In this subsection, we will show the results 
of the above measurements carried on the packings generated by 7 = 1.5 x 10 _3 m/s for 
emulsions and 7 = 2 x 10 4 m/s for granular materials. The measurements done on the 
packings generated by other 7s are basically the same. 
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1. Micromechanical properties 



The normal force distributions in emulsions and granular materials at different pressures 
are plotted in Fig. and respectively as a function of the magnitude of forces normal- 
ized by the average forces. At low pressures, the distributions show a plateau below the 
average force which has been considered as the signature of the jamming transition ("u) ]. 
The distributions show a broad exponential tail for forces larger than the average which ex- 
tend up to eight times the average force for the packings with the lowest pressures. On the 
other hand, V parameters in Fig. |Ht and Fig. have small values at low pressures, which 
means the forces are distributed heterogeneously in space, i.e. they are localized. Combin- 
ing the information from the force distributions and the Y parameter we may conclude that 
at low pressures there are very large forces in the packings and the forces are distributed 



heterogeneously. This agrees with the picture of force chains w 



rich have been visualized in 



granular matter both in experiment |36( and simulations 



13, 



38, 391. These force chains 



sustain most of the external loading. At low pressures, the number of the force chains are 



small and they are well separated 



13j . We also point out that the experiments with confocal 



microscopy of compressed emulsions systems 



20il did not find evidence of force chains in the 



bulk. Our computer simulations of frictionless grains could not reproduce these results as 
we find evidence of localization even for frictionless droplets. 

As the pressure increases, the distributions gets narrow; the exponential tails bend down 
and transform into a Gaussian-like one, as shown by the fitting of the dash lines in Figs |EK 
andEk for grains and emulsions. The deviation from the exponential tail indicates that the 
very large forces (in comparison with the average force) are disappearing. At the same time, 
T increases sharply with pressure in Fig. |Ht and Fig. ^jp, which indicates that the packing 
is becoming homogeneous. This sharp increase happens roughly at 10 MPa for granular 
materials and 50 Pa for emulsions. Similar behavior of T has also been reported in 
for grains and the increase of T was understood as the indication of the disappearance of 
well-separated force chains. 

The distribution at low pressures, containing a plateau and the exponential tail, can be 
fitted by various expressions. Here we choose one from [2^ . 

P(/) = a/ A exp(-(l + A)/). 

This expression comes from a simple Boltzmann equation theory and the power law coeffi- 
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FIG. 9: Distribution of the normal forces in emulsions. Measurements were carried on the packings 
generated by 7 = 1.5 x 10~ 3 m/s. 

cient A is determined by the packing geometry of the system. The fits are shown as solid 
lines in the plots. 

The tangential force distributions in granular materials also show exponential tails and 
their pressure dependence is weaker than in the normal force distributions, as seen in Fig. 
IHb- The distributions of E are plotted in Fig. Eli for five pressures. As the pressure decreases, 
the distribution of S shifts slightly to higher values. This is due to the fact that, in our 
protocol, the system arrests itself more quickly during the construction of high pressure 
packings than low pressure cases. Therefore, the particles in packings at lower pressures 
have more time to slide and have larger tangential displacements; therefore larger tangential 
forces. 

2. Micro structural properties 

The radial distribution function (RDF) describes how the particles are distributed ra- 
dially. In Fig ^2 the radial distribution functions shows a prominent peak at r = ro, the 
second peak at r = 2tq and a sub-peak approximately at r = \/3ro- As the system ap- 
proaches the jamming transition from high pressures, the distribution functions shift to the 
right due to the increase of system size (Figs. [TTk-1 and lllb-1 for emulsions and grains, 
respectively). The first peak decreases in width and increases in hight, as shown in Fig. 
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FIG. 10: r vs pressure in emulsions. Measurements were carried on the packings generated by 
7 = 1.5 x lCT 3 7Tl/s. 

ITTk-2 and [TTb-2. However, we do not find a tendency for the first peak to evolve into a 
5 — function as reported in Q]. In 3^ Silbert et al. reported the power-law behavior 
of Eq. ()15|) near the first peak in granular materials. However, from our data we can not 
clearly identify this power-law region. At low pressures we find r a « 2R, which means that 
the particles barely touch each other. 

The splitting of the second peak in the RDF observed in Fig. [TTb-2 and Fig. ITTb-3 
has been reported by various authors studying dense liquids and granular materials 4lL 
q, y, q The authors of Ref. |4l| consider this splitting to be a structural precursor 
to the freezing transition and the development of long-range order. They also attribute 
this splitting to four-particle hexagonally close-packed arrangements. As shown in Fig. 
the second nearest separation in the hexagonal packings is V3r Q . Therefore the splitting 
sub-peak possibly indicates the appearance of small hexagonal clusters like the one in the 
figure. 

The distributions of contacts are shown in Fig. IT3*k for granular materials and in Fig. IT3*b 
for emulsions. Both distributions can be fitted by Gaussians. As the pressure is increased 
the distributions move to the right and get slightly narrower. The floaters (those particles 
with zero coordination number) exist in granular materials at low pressure (around 12% at 
46.74 kPa) while the fraction of floaters is much lower (2% at 0.498 Pa) in emulsions. As 
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FIG. 11: Radial distribution functions (RDF) in the packings of emulsions and granular materials, 
(a) the complete RDF in granular matter, (b-1) and (b-2) are the first and second peak of RDF 
in granular materials, (c) the complete RDF in emulsions, (d-1) and (d-2) are the first and second 
peak of RDF in emulsions. 




FIG. 12: Illustration of possible origin of the sub-peak. 
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FIG. 13: Distribution of contacts in the packings of granular materials. The lines are Guassian fits 
to the distributions, excluding the floaters. 
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FIG. 14: Distribution of contacts in the packings of emulsions. The lines are Guassian fits to the 
distributions, excluding the floaters. 



expected the number of floaters decreases as pressure increases. Excluding the floaters, we 
can fit the distributions with a Gaussian function, as shown by lines in Fig. IT3k and[*H?b. 
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V. SUMMARY 



We have studied the jamming transition in packings of emulsions and granular materials 
via molecular dynamics simulations. Power-law scaling is found for the vanishing of the 
pressure and excess number of contacts as the system approaches the jamming transition 
from high volume fractions. The emulsion system jams at <p c = 0.645 independent on the 
construction histories while granular materials jam at C between 0.576 and 0.604 depending 
on the construction histories. We found that the preparation protocol has strong effect on 
the tangential forces in granular materials. Longer construction times of the packings allows 
the particles to relax and slide against each other and therefore to have larger tangential 
displacements, which lead to larger tangential forces. Isostaticity is found in the packings 
close to the jamming transition in emulsions and in granular materials in the limit of infinite 
friction and slow compression rate. Heterogeneity of the force distribution increases while 
the system approaches the jamming transition, demonstrated by the exponential tail in the 
force distributions and the small values of the participation number. However, no signatures 
of jamming transitions are observed in structural properties, like the radial distribution 
functions and the distributions of contacts. 
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